function VARbc = doCholVARbootstrap_BiasCorrected_CS(VAR,nboot,clevel )

% Resample full cross-sections using Wild bootstrap. Draw Tx1 vector of
% {-1,1}. Then multiply the residual for all countries in year t by the -1
% or 1 to draw the cross-section.

res      = detrend(VAR.res,'constant');

for jj=1:nboot;
    
        
        % Resample full cross-section of residuals within a time period.
        % Jointly draw all residuals in a given year.        
         Tmax = max(VAR.Ti);                                              
         cal  = (min(VAR.Year):1:max(VAR.Year) )' ;
        rr    = 1-2*(rand(Tmax,1)>0.5);                                    % draw random vector of length Tmax with -1 or 1   
        resb = zeros(VAR.T,VAR.Ny);
        for nn =1:VAR.T;
            resb(nn,:) = rr(find(cal==VAR.Year(nn))) .* res(nn,:) ;        % draw residuals: multiply {-1,1} on residual, same -1 or 1 for all residuals in each cross-section
        end
        
        
        t  = 1;
           Yb = zeros(VAR.T,VAR.Ny);
           Xb = zeros(VAR.T,VAR.Ny*VAR.p); 
           
        for i=1:max(VAR.ii);                                                   
            Xb(t,:) = VAR.X(t,1:VAR.Ny*VAR.p) ;                                
            for j=t:(t+VAR.Ti(i)-1);                                    
                 Yb(j,:)   = Xb(j,:)*VAR.bet_bc(1:VAR.Ny*VAR.p,:) + [VAR.D(j,1:(VAR.N-1)) , 1]*VAR.bet_bc((VAR.Ny*VAR.p+1):end,:)+ resb(j,:);               
                 if j<t+VAR.Ti(i)-1;
                 Xb(j+1,:) = [Yb(j,:), Xb(j,1:((VAR.p-1)*VAR.Ny))];  
                 end
            end 
            t = t + VAR.Ti(i);
        end
        
        % run VAR with this boostrappred data
        VARBS.I_m= VAR.I_m; VARBS.p=VAR.p; VARBS.Ny=VAR.Ny; VARBS.N=VAR.N; VARBS.irhor=VAR.irhor; VARBS.Perm=VAR.Perm;
        VARBS.normalize=VAR.normalize;
        VARBS.Y   = Yb;
        VARBS.X   = [Xb VAR.D(:,1:(VAR.N-1)) ] ;                                 % add fixed effects (dummies) here
        VARBS.bet = [VARBS.X ones(VAR.T,1)] \ VARBS.Y;      
                
        % bias correct estimate
        VARBS.bet   = [VARBS.bet(1:(VAR.Ny*VAR.p),:) - VAR.bc_diff; VARBS.bet(VAR.Ny*VAR.p+1:end,:)];
        VARBS.res   = VARBS.Y - [VARBS.X ones(VAR.T,1)]*VARBS.bet;
        VARBS.Sigma = VARBS.res'*VARBS.res / ( VAR.T - VAR.Ny*VAR.p - VAR.N); 
        % compute Cholesky impulse response function
        VARBS   = irfvar(VARBS);
        
        VARbs.irf(:,:,:,jj)   = VARBS.cholIRF;
             %BS_betas(:,:,jj) = VARBS.bet;
        clear VARBS
end 

% IRF confidence intervals for bias corrected VAR
VARbc.CholirsH1=quantile(VARbs.irf,(1-clevel(1)/100)/2,4);
VARbc.CholirsL1=quantile(VARbs.irf,1-(1-clevel(1)/100)/2,4);
VARbc.CholirsH2=quantile(VARbs.irf,(1-clevel(2)/100)/2,4);
VARbc.CholirsL2=quantile(VARbs.irf,1-(1-clevel(2)/100)/2,4);
VARbc.Cholirsmean=mean(VARbs.irf,4);
VARbc.Xb     =Xb;
end

